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ABSTRACT 

The dynamics of helium shell convection driven by nuclear burning establish the conditions for runaway 
in the sub-Chandrasekhar mass, double detonation model for Type la supernovae, as well as for a 
variety of other explosive phenomena. We explore these convection dynamics for a range of white 
dwarf core and helium shell masses in three dimensions using the low Mach number hydrodynamics 
code Maestro. We present calculations of the bulk properties of this evolution, including time-series 
evolution of global diagnostics, lateral averages of the 3D state, and the global 3D state. We find 
a variety of outcomes including quasi-equilibrium, localized runaway, and convective runaway. Our 
results suggest the double detonation progenitor model is promising and that 3D, dynamic convection 
plays a key role. 

Keywords: convection - hydrodynamics - methods: numerical - nuclear reactions, nucleosynthesis, 
abundances - supernovae: general - white dwarfs 


1. INTRODUCTION 

Despite being rather inert embers of evolved stars on 
their own, white dwarfs (WDs) manage to be at the 
center of many of the Universe’s most spectacular explo¬ 
sions through interactions with companion stars. One of 
the most recently proposed manifestations of such an ex¬ 
plosion are “.la” events (Bildstcn ct al 2007), with the 
decimal point meant to indicate the events are about 
a tenth the brightness for a tenth the time of type la 
supernovae (SNe la). 

The proposed .la progenitor consists of a car¬ 
bon/oxygen (C/O) WD accreting from a helium-rich 
companion. A prominent example of such a system are 
AM Canum Venaticorum (AM CVn) binaries ( 

1995; Nelemans 2005). In their Letter, Bildsten et al. 
(: ') calculate that under the right conditions the ther¬ 

monuclear timescale in an AM CVn’s helium envelope 
can approach the dynamical timescale, possibly estab¬ 
lishing conditions for a detonation which consumes the 
envelope but leaves the WD core intact. This yields 
a relatively faint transient a tenth the brightness of a 
normal SNe la . A unique aspect of these calculations 

1 As is conventional in the literature, we define normal SNe la 
as having lightcurve and spectral properties consistent with the 
dominant population (about 70%, see Li et I (2011)) of observed 
SNe la. Such properties and contrasting peculiar events are dis¬ 
cussed in Branch et al. (1993); Phillips (1993); Hillebrandt et al. 


is the unprecedentedly low ignition pressures, which is 
related to the unprecedentedly low masses of the he¬ 
lium envelopes considered. Previous work considering 
similar systems in the context of double detonations, in 
which the helium detonation triggers a detonation of the 
WD core, assumed higher shell masses (Nornoto 1982a; 
Woosley et al. 1986; Livne 1990; Livne & Glasner 1990, 
1991; Woosley & Weaver 1994; Garcia-Senz et al. 1999), 
excepting a data point in ( 2 ' >) and artificial 

detonations in Livne & \r ( 995). 

As suggested by the authors of the Letter, many took 
on the task of a detailed reexamination of these systems 
with lower mass helium shells. A particularly broad and 
detailed reexamination was carried out by 
Kascn (2011). As they demonstrate, sub-Chandrasekhar 
mass (sub-Mch ) C/O WDs with low-mass helium shells 
can yield a variety of explosive phenomena, including he¬ 
lium novae, double detonations, and deflagrations that 
consumed the envelope, leaving behind a hot core. The 
potential to produce such a variety of transient events 
motivates extensive theoretical inquiry, especially as we 
approach first light for the Large Synoptic Survey Tele¬ 
scope (I zic ' . 2008). A great deal of this inquiry 
has been carried out, with tantalizing results. 

Much of the focus has been on these systems as double 

(2013) 
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detonation SNe la progenitors. Detonation of the C/O 
core appears to be very robustly triggered by compres¬ 
sion waves if detonation occurs in the helium shell ( 
et al. 2007, 2010; Woosley & Kasen 2011; Shen & Bild- 
sten 2014), even in the case of asynchronous, asymmetric 
ignition points ( ). This makes sub- 

Mch promising candidates as SNe la progenitors. Thin 
helium shells have been shown to be capable of carry¬ 
ing sustained detonations, and may even contribute to 
features found in SNe la observations ( 

2012). Synthetic spectra and light curves indicate that if 
the C/O core detonates and dominates over helium shell 
effects in the observables, many sub-Mch progenitor sys¬ 
tems are promising candidates for normal SNe la ( 
et al. 2010; Woosley & Kasen 2011; Fink et al. 2010; 

). In particular, work to date favors 
C/O cores that are near 1.0 M 0 and hot if they are to 
produce normal SNe la. The core-only (no He shell) ex¬ 
plosions of al. (2010) agree best with normal SNe 

la properties for C/O WDs near 1.0 Mq. The ID mod¬ 
els and subsequent synthetic observables of Woosley k 
Kasen (201 ) agree best with normal SNe la’s for their 
“hot” models in which the core relaxed to a luminosity 
of 1.0 L e as compared to their “cool” models, relaxing 
to 0.01 Lq before helium accretion was modeled. 

Further, delay time distribution calculations based on 
binary population synthesis find distributions and rates 
for sub-Mch SNe la progenitor models consistent with 
being at least one plausible dominant channel for re¬ 
producing distributions and rates based on observations 
(Ruiter et al. 201 ). Similar calculations focusing on a 
subset of sub-Mch progenitors find they may be the pro¬ 
genitors of SNe lax (Wang et al. (2013), though see also 
Liu et al. (2015b, a)). Geier et al. (2013) present obser¬ 
vational evidence for both a helium-accreting sub-M^h 
progenitor system and a high velocity helium-rich star 
that matches the expected properties of the unbound 
companion star following a sub-Mch SNe la. Brown 
et al. (2011) analyze a sample of WD binary systems in¬ 
cluding extremely low mass WDs in the context of AM 
CVn binaries and sub-luminous SNe la. They calculate 
merger rates that are comparable to the observed rates 
of sub-luminous SNe la. Drout et al (2013) compare 
observations of SN 2005ek with many possible models, 
including sub-Mch systems. In particular they argue 
that if SN 2005ek did have a sub-Mch progenitor, an 
edge-lit detonation would be the most viable model. In 
the edge-lit scenario, the detonation in the helium layer 
propagates directly into the core, setting off a carbon 
detonation at the core/shell interface. 

We caution that theoretical studies (and the study re¬ 
ported here) have limits and make many assumptions. 
The importance of realistic compositions and convective 
mixing have been made clear (Kromer et al. 2010; Shen 


k Moore 2014; Piro 2015). In addition we note that 
it is currently computationally impossible in any model 
of the full core/shell system to fully resolve ignition of 
core detonation, which occurs on 0.01 - 1 cm scales for 
densities p = 10 7 10 8 g cm -3 . Instead, such work must 
report the critical conditions achieved in a given compu¬ 
tational cell or group of cells and argue the likelihood of 
them achieving ignition of detonation. This challenge is 
in part addressed in Shen & Bildsten (2014), who carry 
out small-scale, fully resolved calculations of detonation 
ignition in regimes relevant to the C/O core of sub-Mch 
systems. They argue that conditions reported in multi¬ 
dimensional studies of the full core/shell system are suf¬ 
ficient for ignition in many cases, though lower mass 
(roughly, below 0.8 Mq) or O/Ne cores are less likely to 
experience ignition. 

A significant uncertainty remains. Only one¬ 
dimensional (ID) models have demonstrated develop¬ 
ment of a detonation in these lower mass helium shells. 
Multi-dimensional work has focused on assuming igni¬ 
tion of detonation and exploring the consequences. In 
this series of papers, we hope to begin to fill this gap by 
modeling the development of ignition and elucidating 
the detailed 3D properties of the system leading up to 
and at the moment of such an ignition. The first paper 
in the series details our methodology, carries out numer¬ 
ical experiments, and demonstrates the development of 
a localized runaway in a model with a 1.0 Mq core and 
0.05 Mq shell (Zingalc et al. (2013), hereafter Paper 
I). In this paper we apply our methodology to a large 
number of models at higher resolution, carry out a new 
numerical experiment, and include many new analyses 
and diagnostics. The purpose of this paper is to 

• expand our methodology to a much larger suite of 
models, 

• explore what broad outcomes and trends we find 
for simple initial models, and to 

• characterize the bulk properties of these models, 
including global 3D structure, ID averages of the 
3D state, and peak global properties such as the 
properties of the hottest cell in the domain. 

This broad exploration lays a foundation for more tar¬ 
geted analysis of a smaller number of more sophisticated 
models. 

We defer a detailed analysis of potential ignition in 
these simple models, including the geometry, timing, 
number, and statistics of localized igniting volumes, to 
the next paper in this series. We also want to be clear 
that one cost of doing such a broad study is that we 
must use simple models motivated by detailed ID mod¬ 
els. Our methodology also fundamentally limits us from 
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modeling the development of a detonation. Instead, we 
are modeling the dynamics we expect immediately prior 
(a few convective turnover times) to any ignition. Thus, 
we model the development of potential seeds of a deto¬ 
nation or deflagration. The ultimate fate of these seeds 
should be the focus of future work. 

To begin, we describe our models and methodology. 

2. METHODOLOGY AND MODELS 

Our simulations are performed using Maestro , a 
finite-volume, adaptive mesh stellar hydrodynamics 
code suitable for flows where the fluid speed is much 
less than the sound speed . Maestro models the flow 
using a low Mach number approximation—sound waves 
are filtered out of the system, but compressibility effects 
due to stratification and local heat release are retained. 
While the code is mature and has been used in state- 
of-the-art astrophysical modeling, it is also under active 
development. al (2010) is the most recent 

and comprehensive in a series of papers describing the 
low Mach number equation set and numerical algorithm. 
The simulations reported here make use of an improved 
low Mach number equation set; this and the associated 
algorithmic changes relative to the model used in 

. (2010) are described in Appendix . The Mae¬ 
stro source code, including all the code needed to run the 
models reported here, and User’s Guide are available in 
the public release. Please see these for the full details of 
the most current form of the algorithm. 

Paper I lays out in detail our numerical strategy for 
modeling sub-Mch WDs with helium shells and demon¬ 
strates the robustness of results to parameter variation. 
Below we review the methodology and describe the con¬ 
figuration of the suite of models used to broadly sample 
the explosive regimes of sub-Mch systems. 

2.1. Microphysics 

We utilize a general, publicly available stellar equa¬ 
tion of state (Timmes & Swesty 2000; Timmes 2008). 

Ions, radiation, degenerate and relativistic electrons, 
and Coulomb corrections are all incorporated. 

Our nuclear reaction network is quite simple for 
the sake of computational efficiency, enabling a broad 
sampling of parameter space (see §2.5). It is im¬ 
portant to note that Wooslcy & Kasen (2011) em¬ 
phasize two crucial reactions for exploring sub-Mch 
systems: 14 N(e", i/) 14 C(a, 7 ) 18 0 (NCO) and 12 C(p, 

2 https://github.com/BoxLib-Codes/MAESTRO 

3 Exactly quantifying how much less is non-trivial. The code 
has been validated against compressible hydrodynamics up to a 
Mach number of 0.2 (Almgren et al. 2006a), and we generally 
abide by a rule of thumb that Maestro is valid up to peak Mach 
numbers of about 0.3. 


7) 13 N(a,p) 16 0 (CagO-bypass). Additionally, Sh n 
& Bildsten (2009) demonstrate the importance of 
14 N(a, 7 ) 18 F (NagF). While we agree these reactions 
are crucial to understanding sub-Mch models, we can 
neglect them for the purposes of exploring the domi¬ 
nant energetics in the pre-ignition burning. For more 
on our rationale for neglecting them and the conditions 
under which they can be neglected, see Appendix 
We employ a simple network consisting of the isotopes 
12 C, 4 He, and 16 0 and the rates for 3 4 He —U 2 C (triple¬ 
alpha) and 12 C(a,7) 16 0 (CagO). CagO is included be¬ 
cause it can allow for the tracing of 16 O production 
which in turn traces the sites of vigorous burning and 
how polluted the shell becomes with burning products. 

Our baseline reaction rates come from 
Fowler (1988), with screening as in Graboske et al. 
(1973); Weaver et al. (1978); Alastuey & Jancovici 
(1978); Itoli e a ( 1). The CagO reaction rate is 

scaled by a factor of 1.7, as recommended in Weaver 
& Woosley (1993); Garnett (1997). Thermodynamic 
derivatives are held constant over a single timestep as 
described in al. (2008). 

2.2. Initial Models 

Maestro evolves both a ID hydrostatic base state and 
a 3D hydrodynamic state. For spherical problems, such 
as the sub-Mch system, this base state is radial. To 
set the initial conditions for our 3D problem we initial¬ 
ize the base state and map that state onto the 3D grid. 
Our sub-Mch initial models are defined by four param¬ 
eters: the mass of the WD core, M\vd, the isentropic 
helium shell’s mass, Mjje, the temperature at the base 
of the helium shell, Tb ase , and the core’s isothermal tem¬ 
perature, T core . At the interface between the core and 
the shell there is a sharp composition and temperature 
gradient following the prescription described in §2.2 of 
Paper I. We generate our own initial models using an 
iterative scheme that enforces hydrostatic equilibrium 
and the values of T core and Tb aS e while converging on 
the given (Mwd, Mn e ). Figure demonstrates a repre¬ 
sentative initial model. 

We expand upon Paper I by adding a new param¬ 
eter test. Most initial models for the simulations re¬ 
ported here use the same transition width parameter S 
as in Paper I: <5 = 50 km. However, we carried out a 
supplemental suite of simulations for model 11030 (see 
§: ) in which all parameters are the same save for a 6 

one-fourth, one-half, and twice the original magnitude 
of 50 km. The quadrupled resolution on a side in this 
paper allows us to resolve sharper transitions than in 
Paper I. Recall that this tanh-smoothing is necessary 
for the problem to be well-posed. A sharp discontinu¬ 
ity in grid-based hydrodynamics makes it impossible to 
demonstrate convergence as it offers no resolvable solu- 


4 


Jacobs et al. 



Figure 1. A representative initial model with Mwd = 1.2 
Mq, M He = 0.05 M @ , T core = 10 7 K, T base = 1.75 x 10 8 K. 
The shaded region is the convection zone. The dashed lines 
from left to right are: the start of the sponge, the anelastic 
cutoff, and the base cutoff density (see §2. !). Top: Temper¬ 
ature (red) and density (blue) profiles. The inset zooms in 
on the sponge and cutoff radii. Middle: Mass fraction pro¬ 
files of carbon (blue) and helium (red). Bottom: Specific 
entropy profile. 


tion to converge to (see Paper I for convergence tests). 
The 50 km value for 6 in the lower resolution models 
of Paper I provided roughly 10 cells of radial resolution 
over which to resolve the transition. The lowest 6 exam¬ 
ined in this paper, 12.5 km, offers similar resolution of 
the transition. In addition, it is not well-known exactly 
how transitions from core to shell are realized in nature. 
Thus, this parameter study has both a numerical and 
physical motivation. 

The top row of Fig. 11 demonstrates the transition for 
11030d0.25 (one-fourth the original 6) and the reference 
model 11030. See § for discussion of results. 

2.3. Grid Structure 

The 3D grid is Cartesian. For all models except one 
(see §2.5) an octant of the sub-Mch WD is modeled, al¬ 
lowing us to capture 3D effects yet achieve much greater 
computational efficiency and explore a large number of 
models. The impact of simulating an octant instead of 
the full star is investigated in Paper I. As we discuss 
in § 1, the higher resolutions and larger model set pre¬ 
sented here introduce complications at the boundaries 
for octant runs with localized runaway. 

The grid is adaptively refined to focus resolution and 
computational power on the regions of greatest interest. 
To study the dynamics of the convection and nuclear 
burning in the helium shell, we refine zones in which 
X He > 0.01 at a density greater than p cutoff (see §2.4). 
To better resolve the shells, in this study we further 
refine cells with temperatures T > 125 MK. We are sat¬ 
isfied with two levels of refinement for models with a 
0.8 Mq core. The mass-radius relationship for WDs 
means these models will have the largest radius and 
consequently the thickest shells spatially. However, for 
models having M core >1.0 Mq the spatial extent of the 
shell is greatly reduced, which is exacerbated by the fact 
that more massive cores have lower mass helium shells. 
Thus, for all such models we add an extra level of re¬ 
finement for a total of four levels (the base grid, which 
we label level one, and three additional, further refined 
levels). 

Figure 2 illustrates the grid we have described. This 
figure outlines grid patches for each of the levels in the 
initial grid for model 10040H (for details about our adap¬ 
tive mesh algorithm and the definition of grid patches, 
see §5 of Nonaka ct al. (201 )). At the coarsest (base) 
level, all octant runs have a 256 3 resolution with a refine¬ 
ment factor of 2 between levels, leading to subsequent 
512 3 , 1024 3 , and 2048 3 effective resolutions within the 
refined patches. We include one full star run (08130F, 
see §2.! ), which has a 512 3 coarse (base) resolution. The 
strong dependence of radius on mass in WDs results in 
a range of physical resolutions Ax ss 2.5 - 15.8 km at 
the finest level. The ID base state’s resolution is not 
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Figure 2. A representative slice of the initial grid for Mwd 
= 1.0 M s , M He = 0.04 M 0 , T c ore = 10 8 K, T baae = 1.85 x 10 8 
K. The different colors indicate grid patches at different levels 
of refinement. Level 1 is the base (coarse) level. 

adaptively refined; instead, it has a fixed resolution of 
five times that of the finest level: 5120 cells. This factor 
of five is first used and discussed in Zingalc ct ■< (2009). 

2.4. Boundaries 

The boundary conditions for our simulations are re¬ 
flecting on the symmetry faces of octant domains (lower 
x, y, and z), and outflow (zero-gradient) on the other 
faces. A full star simulation has outflow boundary con¬ 
ditions on all faces of the domain. 

As can be seen in Fig. , the grid includes a coarsely- 
resolved region well outside the convective zone. This 



serves to keep the convective surface insensitive to 
boundary conditions. A steep drop in density occurs 
at stellar surfaces (as seen in Fig. ). Without modi¬ 
fication, this rapid decline precipitates a rapid spike in 
velocity to conserve momentum. The advantages of a 
low Mach method become negligible if fluid velocities in 
any zone approach the sound speed. Thus much work 
has been put into developing strategies to address steep 
density gradients at stellar surfaces without significantly 
impacting Maestro’s computational or physical validity. 
The details of these treatments can be found in Paper I 
(§2.3), Zingale et al. (2011), and Nonaka et al. (2012). 

Briefly, two density cutoffs are implemented: the 
anelastic cutoff p an eiastic and the low-density cutoff 
Pcutoff (see Fig. 1). For zones with densities below 
Panelastic; Maestro switches to an anelastic-like velocity 
constraint that helps damp velocities (see 
(2 1 )). Density is held constant once it falls to p cu toff, 
halting the steep decline. To prevent impacting validity 
its value is chosen such that the regions with p < p cu toff 
contain an insignificant proportion of the system’s total 
mass. These cutoffs are supplemented with a numeri¬ 
cal sponge that damps surface velocities ( 

108). Cutoff values for all simulations are discussed in 
§2.5. 

This combination of cutoffs, a sponge, and main¬ 
taining a buffer zone in the computational domain be¬ 
tween the stellar surface and the domain’s boundaries 
enables us to study surface convection in detail over long 
timescales without surface effects significantly impacting 
our results. 


Table 1. Model Set 


label 

(M C ore, MHe) 

Tcore 

Tba s e(t = 0 ) 

Phase 

2-max 

AtCfine 

Panelastic 

tfim 

al/(' 7 ”conv) 

outcome 13 


[M®] 

[ K ] 

[x 10 s K ] 

[xlO 5 g cm -3 ] 

[km] 

[km] 

[xlO 5 g cm -3 ] 

min 

avg 

max 


12030 H 

( 1 . 2 , 0 . 03 ) 

o 

00 

1.75 

10.1 

5300 

2.6 

1.29 

2.3 

00 

bo 

12.1 

1 

12030 

( 1 . 2 , 0 . 03 ) 

10 7 

1.75 

10.8 

5100 

2.5 

1.37 

1.5 

6.4 

8.6 

1 

12020 H 

( 1 . 2 , 0 . 02 ) 

O 

00 

1.75 

6.2 

5500 

2.7 

0.80 

3.1 

10.3 

17.1 

1 

12020 

( 1 . 2 , 0 . 02 ) 

10 7 

1.75 

6.8 

5300 

2.6 

0.87 

3.0 

12.4 

19.4 

1 

11030 H 

( 1 . 1 , 0 . 03 ) 

O 

00 

1.85 

5.7 

6600 

3.2 

0.67 

2.7 

10.6 

17.5 

1 

11030 d 0.25 

( 1 . 1 , 0 . 03 ) 

10 7 

1.90 

6.0 

6400 

3.1 

0.68 

12.1 

4.7 

35.9 

1 

11030 d 0.5 

( 1 . 1 , 0 . 03 ) 

10 7 

1.90 

6.0 

6400 

3.1 

0.68 

9.8 

8.5 

29.1 

1 

11030 

( 1 . 1 , 0 . 03 ) 

10 7 

1.90 

6.0 

6400 

3.1 

0.68 

1.9 

26.3 

13.6 

1 

11030 d 2 

( 1 . 1 , 0 . 03 ) 

10 7 

1.90 

6.0 

6400 

3.1 

0.68 

2.5 

32.6 

9.4 

1 

11020 H 

( 1 . 1 , 0 . 02 ) 

10 s 

1.85 

3.6 

6900 

3.4 

0.43 

1.5 

3.8 

9.0 

q 

11020 

( 1 . 1 , 0 . 02 ) 

10 7 

1.85 

3.9 

6600 

3.2 

0.46 

5.2 

19.5 

27.2 

q 

10040 H 

( 1 . 0 , 0 . 04 ) 

o 

00 

1.85 

5.0 

7700 

3.8 

0.58 

7.1 

23.4 

27.6 

q 

10040 

( 1 . 0 , 0 . 04 ) 

10 7 

1.85 

5.3 

7400 

3.6 

0.62 

4.2 

14.9 

18.4 

q 


Table 1 continued 
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Table 1 (continued) 


label 

(M C ore, Mn e ) 

Tcore 

Tbase (t = 0) 

Phase 

■T max 

AtCfine 

Panelastic 

tfir 

al / (^conv ) 

outcome 1 * 


[Mq] 

[K] 

[x 10 s I<] 

[xlO 5 g cm -3 ] 

[km] 

[km] 

[XlO 5 g cm -3 ] 

min 

avg 

max 


10030H 

(1.0, 0.03) 

o 

00 

1.85 

3.5 

7900 

3.9 

0.42 

1.2 

4.6 

7.8 

q 

10030 

(1.0, 0.03) 

10 7 

1.85 

3.8 

7600 

3.7 

0.45 

1.5 

6.5 

8.3 

q 

08130H 

(0.8, 0.13) 

10 s 

1.85 

9.9 

8300 

8.1 

1.15 

0.6 

3.8 

8.4 

l 

08130F 

(0.8, 0.13) 

10 7 

1.85 

10.9 

16200 

15.8 

1.26 

0.2 

1.3 

6.3 

l 

08130 

(0.8, 0.13) 

10 7 

1.85 

10.7 

8100 

7.9 

1.25 

0.6 

4.2 

10.4 

l 

08120H 

(0.8, 0.12) 

10 s 

1.85 

8.8 

8500 

8.3 

1.03 

1.8 

7.7 

10.8 

l 

08120 

(0.8, 0.12) 

10 7 

1.75 

9.6 

8100 

7.9 

1.22 

8.4 

25.7 

27.9 

l 

08050 

(0.8, 0.05) 

10 7 

2.50 

2.6 

10100 

9.9 

0.19 

0.9 

2.7 

7.1 

c 


a see § : for how these values are calculated 

^ outcomes are designated as (1) localized runaway, (q) quasi-equilibrium, or (c) for convective runaway. See text for details. 


2.5. Model Set 

Maestro’s ability to take large timesteps as well as the 
nature of sub-Mch pre-explosive dynamics make a broad 
sampling of the parameter space in 3D computationally 
feasible. What exactly is the parameter space of inter¬ 
est? To determine this we draw on the results of Bildsten 
tl. (2007) and the many studies they inspired. 

The parameters of greatest interest are the core and 
helium shell mass configurations. The motivating ques¬ 
tion is how ignition develops and how it is character¬ 
ized in minimal helium shell mass systems for a range 
of core masses. Figure 2 of Bildsten et al. (200 ) illus¬ 
trates their determination of the minimum shell masses 
for which the nuclear burning timescale is on the or¬ 
der of the dynamical timescale for isothermal cores with 
T core = 3x 10' K. Such a short nuclear burning timescale 
suggests the possibility of thermonuclear runaway even 
for thin helium shells with Mn e < 0.05 - 0.0125M core 
for 1.0 1.2 Mq cores. This work is extended and deep¬ 

ened in subsequent studies, which are largely consistent 
with the essential results of the 2007 work (S 
sten 2009; Brooks et al. 2015). Woosley & Kasen (2011) 
carry out an extensive set of ID sub-M 0 h calculations 
and generate an analogous figure (Figure 19). They in¬ 
clude the impact of varying T core . For M core = 0.7 M 0 , 
runaway can occur with helium shells having ~ 15% of 
the core’s mass, perhaps not sufficiently thin for SNe 
Ia-like spectra. As M core increases to 1.1 M 0 , runaway 
can be achieved with shells ~ 2.25% of the core’s mass 
for hotter cores (T core ~ 7.5 x 10' K), making SNe Ia- 
like spectra more achievable. The bare (no helium shell) 
ID sub-Mch WD detonation calculations of 
(2( )) suggest systems with M core >1.0 M 0 can yield 

observables in reasonable agreement with the range of 
observed normal SNe la while lower M core systems can 


produce characteristics of observed sub-luminous SNe 
la. (a )10)’s 2D calculations also find they can 

produce many characteristics of the range of observed 
SNe la and that core detonation is triggered by shell 
detonations for (M core , Mne) = (0.810 - 1.385, 0.126 
0.0035) M q . 

Given these studies and the uncertainties involved we 
investigate systems with M core = 0.8,1.0,1.1,1.2 Mq 
and a range of shell masses including M^ e = 0.02 — 
0.13 Mq. Our mass configurations are summarized in 
Fig. 3 and compared with the minimum shell masses es¬ 
timated by others. In choosing shell masses we had to 
balance a desire to model low-mass shells near the lower 
limit of models that run away in ID with a need for the 
simulations to be computationally feasible. Lower mass 
cores can take many convective turnover times to reach 
runaway for minimum mass helium shells whereas min¬ 
imum mass shells can be difficult to resolve for higher 
mass cores. As a result, we have the points marked in 
Fig. that track near the ID lower limit but to varying 
extents. 

Due to the importance of T core demonstrated in 
i (2011), we include models with 
Tcore = 10 7 ,10 8 K. Finally, we vary Tb aS e from 175 MK 
to 250 MK. These interface temperatures are intended to 
be roughly what we would expect a few to several con¬ 
vective turnover times before runaway, based on both 
our own numerical experiments and the ID literature. 
Table lists the details of our model set. 

While our models are motivated by the literature they 
are not necessarily likely to be realized in nature and 
are not the result of detailed stellar evolution calcula¬ 
tions. Our focus is on broadly sampling the parameter 
space, characterizing the relationships between param¬ 
eters and possible explosive outcomes, and quantifying 
the salient trends that emerge. This will guide future 
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Figure 3. The crosses are the core-shell mass configurations 
modeled in this paper. For comparison, the shaded region 
is the range of minimum shell masses capable of initiating 
ignition as given in Fig. 2 of ;t< a <:t al. (200 ) (using their 
tnuc = tdyn,10£dyn lines). The hatched region is the range 
of minimum shell masses that yield either a deflagration or 
detonation as given in Fig. 19 of Kasen (2011). 

The lower bound is for their “hot” models, while the upper 
is for “cold” models. 


Ignition: 11030 



work studying particularly interesting parameter config¬ 
urations using more realistic initial models and detailed 
nucleosynthesis. 

3. RESULTS 
3.1. Outcomes 

Our broad sampling of parameter space explores sev¬ 
eral different model configurations. We find a range of 
outcomes for these models: localized runaway, quasi¬ 
equilibrium, and convective runaway. The last column 
of Table denotes the ultimate outcome of each model. 

Localized runaways represent possible seeds of defla¬ 
gration or detonation in the helium shell. All runs in 
this category have localized volumes of fluid that ex¬ 
perience rapid temperature runaway to about 1 GK. 
They also have peak Mach numbers less than 0.3 be¬ 
fore runaway, and less than 0.2 for the majority of the 
simulated time. Fig. plots some of the key properties 
of interest over time for an igniting run (model 11030). 
This plot is representative of the general behavior of 
runs experiencing localized runaway. The plot demon¬ 
strates that the temperature of the hottest cell in the 
domain initially follows a trend similar to that of the 
laterally averaged peak temperature. As the model ap¬ 
proaches runaway the hottest cell increasingly deviates 
from the background conditions. We also find that the 
convectively unstable region moves deeper into the star. 

4 We track the cell with the largest temperature in the entire 
domain, but we caution this is not a Lagrangian measure—it is 
simply the hottest cell without regard to the cell’s mass density 


Figure 4. Temperature, radius, and density over time for 
model 11030. Top: Temperature is plotted in red. The 
solid line is the temperature of the hottest cell in the en¬ 
tire computational domain. The x’s trace the peak later¬ 
ally averaged temperature. The green, cyan, and black all 
plot radii. Green is the radius of the hottest cell in the do¬ 
main. Black is a trace of the core/shell interface based on 
the radius at which the average An e composition is 0.9. The 
cyan plots the base of the convective region and one pressure 
scale height above this base. The inset is a 2D temperature 
slice centered on the site of runaway, demonstrating its lo¬ 
calized nature. The runaway happens at a boundary, hence 
half the inset being white (no temperature data outside the 
boundary). Bottom: The solid line plots the laterally av¬ 
eraged density at the radius of peak average temperature. 
The dashed line is the critical density above which ignition 
is expected according to Eq. . 

This suggests that vigorous burning and the convection 
it drives can result in significant changes in the den¬ 
sity and composition of burning sites (discussed more in 
§3.4). 

As one might expect, the radius of the hottest cell 
moves radially inward along with the base of the con¬ 
vectively unstable region. However, at the end we see 
that this radius moves outward in many models. If con¬ 
vection is able to transport an ignition seed to a signifi¬ 
cant height above the core/shell interface then it makes 
“edge-lit” double detonation models workable. In the 
edge-lit scenario, carbon detonation of the core is trig¬ 
gered by the propagation of the helium detonation wave 
at the surface of the core. This model is generally dis¬ 
favored because it has been shown that it requires the 
initial helium detonation to go off at a substantial height 
above the core/shell interface (Garcla-Senz et al. 1999; 
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Woosley & Kasen 2011). If convection is more effec¬ 
tive than expected at transporting the detonation seed 
then the edge-lit scenario needs to be considered more 
seriously. 

Unfortunately, many of the localized runaway events 
in our models happen near the boundary of the octant 
being simulated. We stress that we have carried out a 
full star run for the localized runaway model 08130 (and 
also in paper I for a 10050 model) and still find local¬ 
ized runaway as well as a radius significantly above the 
interface. We have also carried out simulations in which 
the temperature of cells is limited to be below 3.5 MK 
and see that many localized runaways occur far from the 
boundary, though the initial runaway happens preferen¬ 
tially at the boundary. So while we are confident the 
localized runaway is not a boundary effect, the elevated 
radius of the hottest cell in fig. 4 cannot be ruled out 
as a boundary effect in all runs. This is an important 
issue that will be resolved in the next paper in this se¬ 
ries, which focuses on the timing, thermodynamics, and 
geometry of ignition in this suite of simple models. 

In contrast, quasi-equilibrium simulations balance nu¬ 
clear burning with convective cooling for many convec¬ 
tive turnover times (at least an average of 3.8 or more 
turnovers, see Table 1 for a range of turnover estimates 
and §3.4 for how we calculate these). Fig. 5 demon¬ 
strates this case for a model we ran for a particularly 
long time. We cannot say these runs have reached 
an equilibrium between burning and convective cool¬ 
ing because neither peak nor average base temperatures 
plateau, nor do we model any energy sinks (cooling). If 
we had the computational resources to run these sim¬ 
ulations indefinitely it may be the case they would ex¬ 
perience a runaway. What we demonstrate instead is 
that these models are stable against immediate runaway, 
i.e. runaway within a few convective turnover times. 

We have also included a model similar to 

( 11)’s model 8HB, which experiences a helium 

nova (convective runaway). This model had a higher 
interface temperature than most runs, 2.5 MK instead 
of around 1.9 MK, to facilitate reaching runaway condi¬ 
tions without expending more computational resources 
than necessary. Within the low-Mach limits of Mae¬ 
stro we also find convective runaway, even with the el¬ 
evated interface temperature. As the base temperature 
increases from burning, the turnover rate of the convec¬ 
tive shell is able to increase without plateauing until the 
Mach number of the fluid gets too large for us to track. 
This suggests such a thin shell is able to rapidly trans¬ 
port the energy release of nuclear burning. In contrast 
to localized runaway, this is a more global phenomenon 
and could develop into something like a helium nova, 
as argued in n (20 LI). The time series 

data for this convective runaway is plotted in Fig. ( . The 


Quasi-Equilibrium: 11020 



time [s] 

Figure 5. The same plot as in fig. for model il020, 
demonstrating quasi-equilibrium. 

existence of two regimes, convective and localized run¬ 
away, suggests researchers should investigate the tran¬ 
sition from one to the other. The conditions of this 
transition point will be important for determining the 
minimum helium shell mass capable of achieving local¬ 
ized runaway. 

3.2. Temperature 

While the mass of the core and helium shell play a 
primary role in determining the thermodynamic condi¬ 
tions at the base of the shell, there are secondary de¬ 
terminants. Varying evolutionary histories can result in 
accreting C/O WD primaries of varying temperatures. 
A history of helium flashes may heat the WD surface. 
This enables systems with similar mass configurations 
to have noticeable differences in burning conditions at 
the core/shell interface. 

Our parameterization of the initial model allows us to 
vary the initial temperature of the actively burning base 
of the helium shell. However, in our attempts at varying 
the base temperature we find only a relatively narrow 
range of options can be feasibly explored with our cur¬ 
rent methods. Low initial temperatures will either not 
be able to initiate sufficiently vigorous burning to allow 
a study of ignition or will establish a trend of growing 
average base temperature that will steadily build until 
either ignition or quasi-equilibrium is achieved. How¬ 
ever, the computational resources required to reach a 
dynamically interesting stage of burning with a low ini- 
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Convective Runaway: 08050 



time [s] 


Figure 6. The top plot is the same as in fig. for model 
08050. The bottom plot shows the average velocity magni¬ 
tude in the convective zone as well as the peak Mach number 
in the domain. 

tial base temperature can be substantial. Too high a 
base temperature will either lead to unphysical ignition 
by not allowing time for convection to be established or 
will push the convective velocity beyond Maestro’s abil¬ 
ity to model. Thus, for a given model we choose an 
initial base temperature low enough to allow for convec¬ 
tion to be established and for several convective turnover 
times of evolution, and high enough to reach a scientif¬ 
ically interesting stage of burning while using feasible 
amounts of computational resources. 

The influence of the isothermal core’s temperature 
is easier to investigate with our methods. 

Kasen (2011) find that their synthetic spectra and light 
curves come closest to resembling observations of type 
la’s for their “hot” models in which the accreting C/O 
WD relaxes into thermodynamic equilibrium with a lu¬ 
minosity of 1 L 0 before accretion is modeled. The hotter 
core enables runaway in thinner shells than the colder 
core. This motivates our exploration of our own “hot” 
and “cold” (10 8 and 10' K) models (hot models are in¬ 
dicated with an ’H’ in their label in Tab. 1). 

Our results are consistent with that of 

n (2 111). Hot cores allow for initiation of local¬ 
ized runaway at lower densities for a given core/shell 
mass configuration. This is largely due to an expanded 
core radius in hot runs, and thus a lower density at the 
core/shell interface where the burning occurs. The lower 
density favors higher temperatures at runaway as well. 
In Fig. we compare a hot and cold run to demonstrate 
these phenomena. 


12020H vs 12020 



time [s] 


Figure 7. Comparison of several properties for a hot (10 s 
K isothermal core, solid lines) and cold (10' K, dotted and 
dashed lines) model with a 1.2 M® core and 0.02 Mg shell. 
Top: In red we plot the global peak temperature. In green 
we plot the radius of the temperature peak from a lateral 
average of the 3D data. Bottom: In black we plot the 
laterally averaged density at the radius plotted in green in 
the top panel. 

3.3. Localized Runaway 

We have demonstrated that many of our models 
achieve localized runaway through bulk diagnostics and 
time-series data, comparing them to ID results. This 
answers a major question we are exploring: is ignition 
found in ID codes consistent with 3D models? We argue 
the localized runaway we find is consistent, though do 
caution that localized runaway should not be thought 
of as ignition. The localized runaway reported here may 
ignite deflagrations or detonations, but there is insuf¬ 
ficient evidence and analysis in this paper to make a 
definitive determination. The next paper in this series 
will focus in part on determining the likelihood of such 
ignition. For now, we move to characterizing the local¬ 
ized runaway found in our models. 

The ID studies we have discussed necessarily model 
ignition as simultaneous across a spherical shell. 

( II17) contribute 2D simulations including a vari¬ 
ety of detonation seed geometries, following up later in 
Fink et al. (2010) with 2D simulations seeding a single 
detonation in a larger range of core/shell masses includ¬ 
ing thin shells. Moll & Wooslev (2013) have contributed 
2- and 3D studies in which detonation is seeded at mul¬ 
tiple points with variations in geometry as well as tim¬ 
ing. A key conclusion of these multi-D investigations is 
that detonation in the helium shell very robustly trig¬ 
gers detonation of the C/O core via radially propagating 
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compression waves generated by the helium detonation’s 
shock front traversing the shell . 

Assumptions about how many detonations to seed, 
where to seed them, and their timing impact the ulti¬ 
mate outcome of the double detonation. Single-point 
helium detonations lead to viewing-angle dependences 
not expected in normal SNe la, though the dependence 
becomes weaker for lighter helium shells (Krorncr et al. 
2010). If detonated at a great enough altitude above 
the core/shell interface, helium detonations can directly 
ignite carbon burning instead of detonating indirectly 
with converging shock fronts (Garcia-Senz et al. 1999; 
Wooslcy & Kasen 2011; Moll & Woosley 2013). The 
size and shape of the detonation can also impact how 
easily core detonation can be triggered, especially for 
the thinner shells considered to be the most promising 
candidates for modeling normal SNe la (Moll A Woosley 
2013). 

In light of this, what do our models suggest? To assess 
the conditions that foster ignition we track the hottest 
0.005 % of cells in the computational domain immedi¬ 
ately prior to runaway. Fig. plots histograms of the 
radii and densities of these cells in addition to a spheri¬ 
cal projection of their angular locations for model 08130 
(see Table 1). The spherical projection illustrates the 
two regions in which volumes of hot fluid develop: at 
the base of convective inflows and at the intersection 
of outflows from neighboring convective cells (see §3.4). 
This is determined by contrasting projections like that 
in Fig. with renderings of convective outflow like that 
in Fig. . The density histogram includes a reference for 
a critical runaway density given by Eq. 8 in >oslcy fc 
Kasen (2011): 

Pcr.wK = (1-68 x 10- 4 exp(20/T 8 )) 1/2 ' 3 , (1) 

where T s = T/10 8 is the temperature in units of 10 8 K. 
This is a rough critical density above which violent, hy¬ 
drodynamic runaway is expected and below which con¬ 
vection is expected to be efficient enough to transport 
any energy generated by nuclear burning. For the pur¬ 
poses of this paper, we denote it as a critical density 
above which we expect localized runaway to be possi¬ 
ble. 

This density is based on ID models and thus is only 
fairly compared to lateral averages of quantities in our 
3D simulations. Note that Fig. demonstrates localized 
runaway is achieved even though the hottest cells have 

5 It is important to note the carbon detonation in these studies 
relies upon assuming certain conditions being achieved in a given 
computational cell will lead to detonation. Current studies of the 
full core and shell system do not have the resolution to model the 
initiation of carbon detonation fully self-consistently. See Slicn & 
Bildsten (2014) for detailed carbon detonation calculations. 


densities significantly below this critical density. This 
is not surprising as the cells trace a 3D model. When 
we consider laterally averaged quantities, Eq. is often 
an excellent predictor of ignition. Ignition is achieved 
almost exactly as the average density at the peak burn¬ 
ing region surpasses the (temperature-dependent) criti¬ 
cal density in Fig. , whereas the same average density 
is well below the critical value in the quasi-equilibrium 
case of Fig. . 

This predictor does not work as well for our models 
with 0.8 Mg core masses. These models achieve local¬ 
ized runaway despite being quite a bit below the critical 
density. This could either be an effect of our models 
capturing the 3D dynamics and thus of scientific inter¬ 
est, or it could be a result of our models being toward 
the upper limit of the minimum shell mass the critical 
density is calculated for in our 0.8 Mq models. In addi¬ 
tion, the critical density is only a rough estimate based 
on the outcomes of several ID models. Still, for cores 
> 1.0 Mq it is quite accurate for our models. We will 
further investigate this critical density and determine a 
3D version in the next paper of this series. 

Table lays out key properties at the time of runaway 
for all models that experienced localized runaway. The 
table includes an estimate of the number of convective 
turnover times modeled before ignition as well as the 
lateral averages of temperature and density at the radius 
of peak burning (r pea k)- Note that these are values based 
on the 3D output with a timestamp nearest that of the 
peak temperature. This is why the peak times tend to 
end on round numbers—our 3D state is output at most 
every tenth of a second. 

Table 2. Ignition Conditions 


model ipeak Ipeak/(^conv) Gbase) (Tbase) (Phase) 

[s] [km] [x 10 s K\ [xl0 5 gcm~ 3 ] 


12030H 

118.2 

8.4 

3096.4 

2.180 

13.877 

12030 

120.4 

6.4 

3009.4 

2.178 

14.541 

12020H 

420.0 

10.3 

3261.6 

2.480 

9.424 

12020 

410.0 

12.4 

3163.7 

2.410 

9.963 

11030H 

319.5 

10.6 

3839.8 

2.462 

7.729 

11030 

205.0 

8.5 

3739.1 

2.456 

7.968 

08130H 

161.7 

3.7 

4372.9 

2.139 

11.464 

08130F 

137.4 

1.2 

4251.7 

2.096 

12.218 

08130 

130.6 

3.3 

4259.6 

2.069 

12.106 

08120H 

240.0 

7.7 

4470.0 

2.181 

10.365 

08120 

710.0 

25.7 

4338.7 

2.041 

11.221 


3.4. Convection 
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Figure 8. Hotspot properties for model 08130 at t = 130.0s. Upper left: histogram of the radii of the hotspots, with bin sizes 
chosen to match the finest level of 3D resolution. Each bin is color-coded to indicate the maximum temperature in the bin. The 
location where helium’s laterally averaged mass fraction is 0.9 and the base of the convective envelope are indicated. Upper 
right: histogram of the densities of the hotspots. Again, bins are color-coded based on maximum temperature. The value of 
Eqn. is indicated. Bottom: projection of the hotspots’ angular location onto a sphere corresponding to the average radius 
of the hotspots. Hotspot pixel sizes correspond to the finest level’s physical resolution, and where hotspots overlap preference 
is given to the one with the highest temperature. The extents of the temperature color bar’s bins are set such that each bin 
contains an equal number of hotspots. 
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The ignition characterized in the previous section is 
fostered by the complex interplay between nuclear burn¬ 
ing and convective dynamics. A detailed understanding 
of convective evolution is crucial, as illustrated by the 
conflicting results of Fink et al. (2010) and Woosley & 
Kasen (2011). 

Fink et al. (2010) get their initial ID models from 
Bildsten et al. (2007), who assume a fully convective 
shell all the way up to the point of ignition. Woosley 
& Kasen (2011) employ a time-dependent convective 
model based on mixing-length theory, allowing for con¬ 
vection to “freeze out.” If the runaway timescale for a 
volume of fluid at the base of a convective zone, in the 
absence of cooling, becomes smaller than the convective 
turnover timescale, convection is no longer able to cool 
the helium-burning layer with the same efficiency. It 
starts to “freeze out.” A fully convective shell requires 
a larger temperature at its base to achieve runaway than 
with one in which freeze-out is allowed. The tempera¬ 
ture differences translate into entropy differences. In 
turn, the lower entropy of the cooler base makes possi¬ 
ble larger base densities in the Woosley & Kasen (2011) 
models than those of ink et al (2010), even when they 
are modeling similar core and shell masses. These den¬ 
sity discrepancies lead to significant discrepancies in ex¬ 
plosive burning and the products yielded. 

Our analysis of convective dynamics begins by estab¬ 
lishing the broad context. Fig. plots volume render¬ 
ings of radial velocity for several models. These act 
as a proxy for convective plumes. Higher core masses 
result in more compact systems with smaller pressure 
scale heights. Thus we see the typical size of a con¬ 
vective cell grows inversely proportionally to core mass. 
The hotspots plotted in Fig. occur at the base of con¬ 
vective plumes of cool in-falling fluid and at regions in 
which outflows of adjacent convective cells collide. The 
in-falling matter increases the density, setting off more 
vigorous nuclear burning. A competition between the 
nuclear burning rate of a volume of fluid and the rate 
at which it can be transported and cooled by convective 
flow is established. 

In our models we see convective dynamics have two 
key impacts. First, convective overshoot serves to push 
the region of active burning deeper into the star, thus in¬ 
creasing the ambient density of nuclear burning sites and 
altering the ambient composition. Second, convection’s 
ability to respond to increasing nuclear energy genera¬ 
tion breaks down for sufficiently energetic burning. This 
is the freeze-out discussed in Woosley & Kasen (2011). 

Convective overshoot is demonstrated in the top plot 
of Fig. 10 for model 11030. The first process to hap¬ 
pen in a model is for convection to be established in 
response to nuclear burning, the bulk temperature pro¬ 
file, and our initial velocity perturbations (see Fig. 1). 


This is seen in the plot of the convective timescale at 
the bottom of Fig. . Initially it experiences a steep 
drop off until stabilizing as convective cells are estab¬ 
lished. After this we see that the location of the convec¬ 
tive base steadily moves radially inward. In response, 
the radius of peak burning moves deeper into the star 
resulting in increased ambient density as well as changes 
in ambient composition. We find in most models that 
experience localized runaway, peak burning radii tend 
to stabilize near the location where composition is 90% 
helium (shell material), 10% carbon (core material). In 
addition, models in quasi-equilibrium do not have such 
substantial radial migration of the convective base. 

Freeze-out is demonstrated in the bottom of Fig. 10. 
Here we plot two different timescales. To calculate a 
minimum nuclear burning timescale we invert the burn¬ 
ing rate, dX i /dt 1 for the most rapidly burning species i 
(for these models, helium). This is done for all radius 
bins r in our lateral averaging. The minimum value of 
this radial slice is used as our nuclear timescale r nuc . In 
sum, 



The other timescale is a conservative estimate of the 
convective turnover time. Again using laterally aver¬ 
aged data, we invert the velocity magnitude (|U|) and 
integrate over the convective region, 


^"conv 



(|U(r)|) 1 dr. 


( 3 ) 


rb is the smallest radius at which the radial entropy pro¬ 
file satisfies ds/dr = 0 and r t is largest radius satisfying 
this condition. The region is shaded in Fig. . This av¬ 
erage value is reported in Table . In addition, Table 
includes a calculation of the minimum and maximum es¬ 
timates for turnover times by simply dividing the length- 
scale r t — rb by the maximum and minimum velocities in 
the convecting region, respectively. We want to stress 
these are conservative estimates. Not all plumes that 
form will extend the full lengthscale, and there may be 
plumes with faster fluid than in other plumes. 

In the timescale plot of Fig. 10 we see that as 
model 11030 approaches runaway the nuclear burn¬ 
ing timescale drops more rapidly than the convective 
turnover timescale. We note that the focus is not on the 
magnitudes of the timescales but their changes. The 
highly nonlinear nature of nuclear burning make com¬ 
parisons of the rate’s magnitude at any given time with 
other timescales uninformative. The changes however 
suggest the nuclear burning rate is shrinking faster than 
the convective turnover rate as runaway is approached, 
suggesting a degree of freeze-out. 

In contrast, quasi-equilibrium models and the convec- 
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Figure 9. Volume rendering of radial velocities for some representative models. Maximum velocities are 2.8, 1.3, 2.0, and 1.7 
xlO' cm s -1 for models 08120H, 10040H, 11030H, and 12020H, respectively. 


tive runaway model have relatively flat peak burning 
radii. Their convection is able to respond to nuclear 
burning before any localized runaway can develop. 

3.5. Varying S and Comparison with a ID Model 

As discussed in §2.2, we have carried out an additional 
numerical experiment exploring the impact of varying 
the 5 parameter that determines the sharpness of the 
transition from core to shell (see Paper I for details). 
In Fig. we plot the temperature profiles of models 
11030d0.25, 11030d0.5, 11030, and 11030d2, which have 
identical initial models except for their 5 parameters: 
12.5, 25, 50, and 100 km, respectively. In addition, 
we plot the temperature profile of model 10HC from 
Woosley ( 111) at a time near ignition (model 

data courtesy of Woosley, private communication). The 


radius of 10HC is offset by 500 km to facilitate com¬ 
parison of transition widths. We plot both the initial 
temperature profile and those from later times. Initially, 
our default 6 value results in a more gradual transition 
than in 10HC. At later times, the temperature profile 
develops a more pronounced peak like that of 10HC, 
though the transition continues to be less sharp than 
10HC and leaves a substantial amount of hot fluid be¬ 
low the temperature peak. We see that our smaller 6 
values, in particular that of 11030d0.25, are about as 
sharp as that of 10HC. 

One potential worry with the thicker transition is that 
it allows a thin shell of hot fluid (80-90% of the peak in¬ 
terface temperature) to exist below the region initially 
unstable to convection (see Fig.I). If convection is not 
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Model 11030 Convection Timeseries 




Figure 10. Top: Here we plot: the radius of the cell with 
the largest nuclear energy generation in the entire domain, 
r(e P eak), in cyan; the radius of this same quantity for the 
laterally averaged data, also in cyan with “x” markers indi¬ 
cating the timestamp the data was calculated for; the radius 
at which the lateral average of helium’s mass fraction (A'He) 
is 0.9, in magenta; the radius of the convective base as deter¬ 
mined by the radius at which the entropy flattens out (see 
Fig. 1), in yellow; and the radius at which the lateral average 
of temperature (T) peaks, in black. Bottom: A comparison 
of a nuclear burning timescale (magenta) with a convective 
turnover timescale (cyan). See text for details. 

well-established in these hot shells as the model ap¬ 
proaches runaway then it becomes hard to determine 
if the runaway is a result of the uncooled hot shell in 
some initial models or the convective dynamics. To ad¬ 
dress this concern, we plot slices of temperature at t = 0 
and late times for 11030d0.25 and 11030 in Fig. . This 
plot includes white vertical lines marking the radius at 
which the laterally averaged velocity magnitude is 25% 
and 50% of its peak value. This gives the reader an 
idea of how much convective cooling penetrates. The 
late-time temperature profile develops into a thin, hot 
layer for both S values, with similar velocity penetration. 
This demonstrates that even with a thicker transition 
the cooling in our models penetrates to the thin shell of 
vigorous burning. The primary impact of the smaller 6 
is to shift the radius at which the thin, hot layer devel¬ 
ops. This is an expected imprint of the initial model, as 
the thinner S locates a more concentrated shell of hot 
fluid at a lower radius, as seen in Fig. . In addition, 
the thicker S results in a larger region of intermediate 
temperature below the thin, hot layer. Finally, we note 
that a history of helium flashes can heat the base of the 


convecting envelope. Nature may in fact realize config¬ 
urations in which hot fluid exists below the convectively 
unstable region. 

We find that all of the runs with varied J’s also ex¬ 
perience localized runaway. In addition, the convective 
dynamics reported in previous sections are qualitatively 
insensitive to varied <5. Figs. and illustrate this us¬ 
ing model 11030d0.25. We also see a new result: the dra¬ 
matic rise of the typical radius at which helium’s mass 
fraction satisfies X^ e = 0.9. This suggests substantial 
mixing as carbon is dredged up, displacing helium. This 
will be explored in more detail in the next paper in this 
series. 

4. CONCLUSIONS AND DISCUSSION 

We have explored 18 physically motivated, simple 
models of convective burning in sub-Chandrasekhar car¬ 
bon/oxygen white dwarfs with an accreted, low-mass he¬ 
lium shell, as well as 3 supplemental models exploring 
a new numerical experiment of the impact of an initial 
model parameter. These systems have been modeled 
in 3D with detailed microphysics using the low-Mach 
hydrodynamics code Maestro. These are the first 3D 
simulations of this phase of evolution for such a broad 
suite of models, and we have drastically expanded the 
diagnostics and analysis in comparison to Paper I. 

As we make clear in §2, our models are simple and 
limited in the number of species and reactions we track. 
This is sufficient for our focus on exploring a large num¬ 
ber of model systems and capturing the dominant ener¬ 
getics and general dynamical trends. Our findings here 
will serve as the foundation for future work focused on 
a smaller number of more detailed models. From this 
set of simple models we can draw several important in¬ 
sights into potential double detonation SNe la or .la 
progenitors as well as other runaway events involving 
sub-Chandrasekhar white dwarfs with envelopes of ac¬ 
creted helium, such as possible helium novae. 

We find that localized runaway is indeed achievable 
in 3D. Whether this develops into a deflagration or det¬ 
onation will be the focus of future work. In particu¬ 
lar, it would be valuable to make use of the Maestro- 
to-Castro mapping developed in Malone (2014). 

Like Woosley & Kasen (2011), we find that hotter cores 
allow for localized runaway to develop at lower densities 
than in models with cooler cores. Thus, hotter cores 
are able to achieve runaway with thinner (lower mass) 
helium shells than cooler cores of the same mass. 

Our results indicate that the complex dynamics of 3D 
convection should not be neglected in investigations of 
double detonation progenitor models. In cases of local¬ 
ized runaway, we find evidence of convective overshoot 
and a steady freeze-out of convection. These effects sub¬ 
stantially impact the density, temperature, and compo- 
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(c) 11030d0.25 t = 375 s (d) 11030 t = 180 s 

Figure 11. Zoomed-in X-Y slices of temperature at the base of the convecting helium envelope for models 11030 and 11030d0.25 
at initial and late times (see subplot labels). Cell edges are overlayed. The two white vertical lines mark, left-to-right, the radius 
at which the lateral average of the velocity magnitude is 25% and 50% of its peak. 


sition of the region in which localized runaway is initi¬ 
ated. 

Future papers in this series will characterize the geom¬ 
etry, thermodynamics, statistics, and timing of ignition 
in greater detail, carry out more realistic simulations 
making use of ID stellar evolution for initial models and 
larger reaction networks, and map our low-Mach results 
into the compressible code Castro to track the develop¬ 
ment of deflagration or detonation. 
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Initial Model 




Figure 12. Comparison of temperature profiles for varying 
5. Top: The temperature profiles of initial models with 
varying 5 along with a realistic ID reference (model 10HC 
from Wooslev & Kasen (2011), radius shifted —500 km to 
facilitate comparison). Bottom: The same reference model 
along with profiles for the laterally averaged temperature 
data from models 11030d0.25, 11030d0.5, 11030, 11030d2 
after ~ 150s of evolution. 
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Figure 14. Same as Fig. for model 11030d0.25 
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Figure 13. Same as Fig. for model 11030d0.25 
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APPENDIX 

A. AN IMPROVED LOW MACH NUMBER EQUATION SET 

In this Appendix, we discuss the improved low Mach number equation set used for the simulations presented in this 
paper, and highlight the differences from the algorithm presented in Nonaka et al. (2010). We also discuss the effects 
of this improved model on some standard test problems. 

Recently, both (.. 2) and . (21 ) introduced a new term in the vertical momentum 

equation for low Mach number stratified flows that enforces conservation of total energy in the low Mach number 
system in the absence of external heating or viscous terms. To compare the new system with the origianl Maestro 
equations, we first define p*,U*,p*, and T*, as the density, velocity, pressure, and temperature, respectively, in the 
low Mach number system, in keeping with the notation used in ( >89) and (. ). The 

perturbational quantities, p 1 '* = p* — po and p 1 ’* = p* — po, are analogous to those in the compressible system. The 
quantities defining the starred fluid approximate, but are not identical to, the quantities defining the fully compressible 
fluid. 

The following equations were derived in (:. )6a) for low Mach number stratified flow with a general 

equation of state: 


with the constraint of the form, 


dp* 

Ik 


+ V • (p*U*) = 0 


sir 

~df 


+ u* • vu* + —v P '’* 

p* 



V • (A)(r) U*) = f3 0 aH , 


(Al) 

(A2) 


where 

A(r) = *“p(i FtUP)*') 

and er = PT/{pc v Pp) ( a 2006b), which in the case of an ideal gas reduces to l/(c p T). Also, we define Pi 0 

as the lateral average of the first adiabatic exponent, 


Pi = d{\ogp)/d{logp)\ s 


(A3) 


( ) modify the vertical momentum equation by noting that a low Mach number representation 

of total energy is conserved if, in moving from the compressible to the low Mach number system, one substitutes 
1/p —► 1/p* — (l/p*) 2 (dp*/dpo)\ s (p* — po) instead of 1/p —► 1/p* in the momentum equation, resulting in 


<9U* 

dt 


+ U* • VU* + — S/p’* = — 

p* p* 


dp* 

dpo 



(A4) 


( ) instead construct the additional term by deriving the low Mach number momentum equation from 

Lagrangian analysis, starting from conservation of total energy. The momentum equation given in ( 013) 

uses the assumption from lmgren cl (2006a) that Tj can be approximated by Pi 0 , and can be written in the form 


<9U* 

~dT 


+ u* • VU* + 




* 9 e r 

p 


(A5) 


We note that ( ) is analytically equivalent to ( ) when Ti = r lo . 

An essential component of the solution procedure for low Mach number equation sets is solving a variable coefficient 
Poisson equation for the perturbational pressure. The Poisson equation can be derived by substituting the constraint 
into the divergence of the momentum equation, 

y - (^ v y’*)= v -(/9°A*)-/3°^^, (A6) 


with 

(p* - Po) 

_ § 


A* = -[U* • VU*] 


P 
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Using their formulation, ( ) show that one can instead solve 


V • 



= V ■ 


(PoA*) - fa 


d(aTL) 

dt 


(A7) 


for p'’*. Analytically, this looks very similar to the original equations, and numerically, it allows one to reuse the 
original solver, simply modifying the coefficients (P%/p* as opposed to P 0 /p*) and the interpretation of the variable 
being solved for (p'’*/Po as opposed to p’’*). While exactly the same solver can be used to solve Equation ( ) as to 

solve Equation ( >), we note that solving for p 1 '* in the new constraint may take more computational effort. This 

arises because in Equation (. 6) the coefficients of p 1 ’* are close to one, since Po is a density-like variable that is close 
to po (and in fact identically equals po for an isentropically stratified atmosphere). The coefficients in Equation (. 7) 
are, by contrast, similar to po (r), which can have large variation over the scale heights within a single calculation. 

We note also that the simplification allowed by writing the equation for p' 1 * in the form Equation (i 7) as opposed to 
more elaborate formulations used in ( ) occur under the assumptions that the flow is non-adiabatic 

and the base state is constant in time. The implications of adiabaticity and a time-varying base state are subjects of 
future work. 


A.l. Test Problems 

Here we look at the impact of the different formulations of the momentum equation on two test problems we have 
used in the development of Maestro as well as one of our science cases. In (:2 ) a study of linear gravity 

waves using the two different formulations as implemented in the Maestro code demonstrated that the eigenfunctions, 
energy conservation properties and pseudo-energy conservation were as expected. We refer the reader interested in 
those results to Vasil et a (2013). 

A. 1.1. Reacting Bubble Rise 

In ( ) we showed an example of three burning bubbles in a stratified white dwarf atmosphere 

(modeled in a plane-parallel geometry, for convenience). The bubbles are buoyant and rise and roll up as nuclear 
reactions release heat in the hotter locations. Here we revisit that problem with the modified momentum equation. 
Figure shows a comparison of the temperature field (T* = T(po, e*)). for the simulation using the new formulation 
vs. the original equations. As we can see, the new formulation generates slightly wider bubbles. 

This difference in width is consistent with an observation made in (2006a), where we compared low 

Mach number simulations of non-reacting buoyant bubbles with compressible simulations. There we saw variation 
between the bubbles as calculated with the different methods (including between the two compressible formulations), 
but observed that the bubbles that evolved via the low Mach number equation set were consistently narrower than the 
fully compressible bubbles. The fact that the new formulation generates slightly wider bubbles suggests in a general 
sense that the solution with the new formulation is closer to the frilly compressible solution than that generated with 
the original formulation. 

This test problem is distributed with the public version of Maestro as reacting_bubble. 

A. 1.2. Internally Driven Convection 

A simple test of convection driven by an analytic heating source was shown in ( )( )8) and again in 

Nonaka et al. (2010). The details of the setup and heating term are given in Section 4.3 of (2008). 

The basic idea of this problem is that there is a convectively unstable region bounded above and below by stable 
regions. An analytic heat source at the base of the unstable layer drives convection. The top of the convective region 
essentially spans to the top of the atmosphere, up to the point where the density drops off steeply. We considered the 
convectively unstable layer to be the “region of interest” for comparison with results from a compressible code, and we 
found good agreement between the original Maestro equations and the fully compressible solution. Here we compare 
results from simulations using the new formulation against both the original results and the compressible results. The 
compressible comparison for this paper is done with the CASTRO code ( 2010). 

For the calculations presented here, as in the previous two works, a plane-parallel approximation is used for the stellar 
atmosphere, and 320 x 512 cells are used to cover the the domain spanning 2.5 x 10 8 cm in the horizontal direction 
and 4.0 x 10 8 cm in the vertical. In Figure we show plots of T*, p 1 ’* /po, and Mach number for the the original and 
new algorithms. Two things are immediately apparent. First, in the convective region, roughly between 10 8 cm and 
2.5 x 10 8 cm, we see good agreement between the original algorithm and the new version. The temperature field and 
Mach number in that region agree nicely. However we also note that above the convective region, the solutions differ. 
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Figure Al. Comparison of burning bubbles with the new formulation and the original formulation of the Maestro equation 
set. 


This is precisely the location of the surface of the star, so the density drops sharply to the cutoff value here. We notice 
that in the new model, the Mach number is much lower in this surface region and p' /po is also much lower (for the 
latter, we note that the plots use different scales to bring out the details in each case). While these differences do not 
appear to influence the behavior of the convective region, the large M in this surface region does have implications for 
the timestep. Overall, the new model seems better behaved in the surface layer. Looking at the stable region beneath 
the convective layer, we do not see much of a difference between the two models. 

As we did when we previously looked at this test problem, we present two diagnostics: the lateral average of the 
temperature, T*, 


( T )j 


1 

iv* 


-l" x 

Ex, 


2=1 


(A8) 


and the RMS fluctuations: 


m, = 




Nr 



(A9) 


where N x is the number of cells in the lateral direction. Figure shows the profiles of (T) and ST /( T) as a function 
of height. In the region of interest, the convective layer, we see that all the solutions agree strongly. The differences 
outside of the convective layer also seem minimal. The increase in T that is seen above the convective layer for the 
compressible solution is because in the compressible code, the material above the atmosphere is not in hydrostatic 
equilibrium (it has reached our cutoff density) and rains down onto the surface. This does not have a large dynamic 
effect because there is not much mass there. 

This test problem is distributed with the public version of Maestro as test_convect. 
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t = 5 5 


Figure A2. T*, p'’*/po, and Mach number for the convection test problem using the original (top) and new (bottom) algorithms. 
Note due to the large difference in p’’*/po between the two runs, the original formulation is clipped at the upper value of the 
plot range. In reality, p'’*/po is two orders-of-magnitude larger in the original formulation than in the new energy formulation 
in the region above the atmosphere. 
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Figure A3. (T) and ST / (T) as a function of height for the convection problem, comparing the original Maestro implementation, 
the new energy formulation, and the compressible solution from CASTRO. 

B. NUCLEAR REACTION ENERGETICS 

To demonstrate that we utilize the smallest set of reactions and isotopes capable of yielding the dominant energetics 
we calculate the rates for relevant reactions using MESA’s 1 flexible and extensive reaction network module (I ton 
et al. 2011, 2013). MESA uses multiple sources for reaction rates, including (20( 6); Angulo < . 1 (' 99); 

(1 8). As we do in our own Maestro network, we scale the 12 C(a,7) 16 0 (CagO) reaction by 

1.7 ( ver Sz Woosley 1993; Garnett ). The code, data, and configuration for these MESA calculations are 
distributed in the Maestro repository. 

To explore how the reactions highlighted by others (see §2. ) might contribute to the energetics leading into ther¬ 
monuclear runaway, our MESA calculations use a composition that allows all the relevant reactions to occur and is a 
reasonable reflection of what we expect in nature. We assume a core of 49.5% 12 C and 16 0, 1% 22 Ne, and a shell of 
99% 4 He, 1% 14 N. Since the energy release of 14 N(e _ , u) 14 C(a, 7) 18 0 (NCO) is dominated by the alpha capture on 
the intermediate 14 C, 10% of the 14 N is converted into 14 C. The site of active burning is assumed to be 10% core 
material, 90% shell material as it is in the simulations being reported here. Finally, all mass fractions are reduced 
by about 5% to allow for 5% of 4 H so that the 12 C(p, 7) 13 N(a,p) 16 0 (CagO-bypass) reaction can proceed. The 
resulting mass fractions are Xi H = 0.05, X4 He = 0.84645, Xi 2 C = 0.047, Xi4 N = 0.0077, Xi4 C = 0.00085, Xisq = 0.047, 
X 22 Ne = 0.001. In the pre-explosive dynamics leading up to thermonuclear runaway the sites of nuclear burning have 
T ss 200 MK, p ss 10 5 10 6 g cm~ 3 (see §2.2). The most energetic reaction rates for these conditions and the given 


http://mesa.sourceforge.net/, version 7503 
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T = 200 MK 



P [g cm 3 ] 

Figure B4. Energy generation of dominant reactions 

composition are plotted in Fig 

The forward NCO reaction does not become faster than the back-reaction until about pa 4 x 10 6 g cm~ 3 and 
even then is only marginally faster for the densities considered here, so it is not included in the figure. A much more 
important consequence of a shell with 14 N is the quite energetic NagF reaction which can generate an appreciable 
fraction of the amount of energy generated by the triple-alpha reaction. The triple-alpha reaction still dominates and 
the dynamics we demonstrate here will be similar with or without the NagF reaction, except that ignition may be 
achieved more quickly with it. For simplicity, we neglect it in this study. A future study building on this work will 
explore more detailed nucleosynthesis and the impact of using a larger reaction network. 

The CagO-bypass reaction is incredibly energetic. We can only justify neglecting this reaction in our calculations 
because we are considering pure helium accretion with no free protons before thermonuclear runaway conditions of 
T al GK are achieved. As discussed in i & Bildste i (2009), even in the event of no protons in the envelope, once 
conditions of T a 1 GK are met alpha-chain reactions of 24 Mg and higher can yield intermediate protons, and the 
14 N(a, 7 ) 18 F(a, p) 21 Ne reaction directly yields protons (in the pre-explosive regime the 1S F is unable to alpha-capture 
at an appreciable rate). Thus the CagO-bypass reaction must be included for any studies of the sub-Mch ignition 
regime. 

From these rate calculations we conclude the triple-alpha reaction dominates energy generation in the pre-ignition 
conditions and expected composition of C/O WD helium-accretors. 
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